function discard  = fsb_det_artifact_mod(conc_idat,filter_param,what)

% FSB : Function to detect and label defective volumes
%
% EXAMPLE:
% discard  = fsb_det_artifact_mod(conc_idat,50,1);
%
% INPUT
% conc_idat:    4D Image data
% filter_param: filter parameter for artefacts.
%               (number between 1 and 50 that determines the
%               amount of noise that the algorithm tolerates).
%               1 is a very strict filter, 50 lets everything pass
% what:         if additional information should be used for the filter,
%               set this to 1
%
% OUTPUT
% discard:     vector of volumes to discard
%
% CALLED BY:
% FSB.m
% fsb_meancalc.m
%
% NOTES
% integral component of fMRI Sandbox. Especially useful for segmented data
%
% Copyright Steffen Stoewer 28/11/07
%
% Copyright 2010 MPI for Biological Cybernetics
% Author: Steffen Stoewer
% License:GNU GPL, no express or implied warranties
% 
% $ Revision 1.0
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

if nargin<3
    what = 0;
end

[a,b,c,d]=size(conc_idat);% get image dimensions

slicex = int16(round(a/2));
slicey = int16(round(b/2));
slicez = int16(round(c/2));

correction = 0.1;

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Define areas inside and outside of brain
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
brain = conc_idat(slicex-round(slicex/4):slicex+round(slicex/4),...
    slicey-round(slicey/3):slicey+round(slicey/3),...
    slicez-round(slicez/2):slicez+round(slicez/2),:); % Takes a couple of slices in the middle from front to back

roi_out1 = conc_idat(:,slicey+round(slicey/1.5):b,...
    1:slicez+round(slicez/3),:);%outside brain
roi_out2 = conc_idat(:,1:round(slicey/4),1:slicez+round(slicez/3),:);%outside brain
roi_out = [roi_out1 roi_out2];

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Determine means of in and out of brain areas
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

brain_avg = fsb_meancalc(brain,4);
roi_avg = fsb_meancalc(roi_out,4);

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Calculate ratio and threshold
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

verhalt = roi_avg./brain_avg;
filter = fsb_avg(verhalt)+(filter_param/100)-correction;
discard_filt = find (verhalt>=filter);

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Check if more stringent threshold desired
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

if what == 1
    diff_idat = avg(repmat(conc_idat(:,:,:,1),[ 1 1 1 d]) - conc_idat)*10;
    discard_diff = find(diff_idat<mean(diff_idat));
    discard = union (discard_filt, discard_diff);
else
    discard = discard_filt;
end

end
